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To distinguish between regular and cha otic orbits in Hamiltonian systems, the Global Symplectic 
Integrator (GSI) has been introduced Libert et al. . 201dj | . based on the symplectic integration 



of both Hamiltonian equations of motion and variational equations. In the present contribution, 
we show how to compute efficiently the MEGNO indicator jointly with the GSI. Moreover, we 
discuss the choice of symplectic integrator, in fact we point out that a particular attention has to 
be paid to the structure of the Hamiltonian system associated to the variational equations. The 
performances of our method is illustrated through the study of the Arnold diffusion problem. 

Keywords: Symplectic integration, chaos, variational equations, MEGNO, Arnold diffusion. 



1. Introduction 

Several detection techniques exist in order to study the regular or chaotic behavior of orbits of Hamiltonian 
systems. All Lyapunov-like methods are based on the resolution of the so-called variational equations giving 



the e volution of deviation vector s assoc iated to a given orbit. These are es sentially the FLI IFroeschle et al.. . 
1997tl. MEGNO ICincotta et all 12003] . SALI [Skokosl . l200ll ] and GALI [Skokos et all [20071] methods. In 
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Libert et all l2010j ]. we have introduced the Global Symplectic Integrator (for short the GSI), a method 



allowing to solve both Hamiltonian equations of motion and variational equations using a totally symplectic 
integration scheme. Based on a comparison of the GSI with non-symplectic integration schemes, it was 
clearly shown that the GSI was more accurate in the detection of regular and chaotic orbits, using the 
SALI chaos indicator, and less time-consuming than non-symplectic methods. Such a kind of symplectic 
integration scheme for the variational equations has also been recently and independently proposed in 
[Skokos et alLl201Cj |. 

The purpose of this work is to show that other chaos detection techniques, for instance the MEGNO, 
can be used jointly with the GSI and their efficiency thus improved. The computation of the MEGNO 
needs to study the time evolution of a single deviation vector. This advantage, with respect to the SALI 
that needs the knowledge of the evolution of two deviation vectors, becomes critical when considering 
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a large amount of orbits on long time spans. Moreover, the use of the GSI does not rely on the use of a 
specific symplectic integrator scheme. Hence, our second aim is to discuss and compare different symplectic 
integrators. It turns out that methods adapted to the structure of the Hamiltonian equations of motion 
are not necessarily suited for the associated variational equations. 

While the study of the dynamics of two well-k nown Hamiltonian systems, Henon-Heiles and the re- 
stricted three-body problem, have been addressed in [Libert et all 120101 ] . we hereby propose to use the GSI 



to improve the detection of slow diff usion in Hamilton ian systems, the so-called Arnold diffusion, analyzed 
according to the model proposed in [Lega et all 1200.^ . 



The organization of the paper is as follows. In Section we briefly summarize the GSI method, and we 
review two classes of symplectic integrators. Section [3] is devoted to the definition of the MEGNO indicator 
and the algorithm used to compute it. An application to the Arnold diffusion problem is then presented in 
Section HI for which the the GSI method is compared to a non-symplectic scheme. Finally, in Section [5] we 
sum up and draw our conclusions. 



2. Global Symplectic Integrator 
2.1. Method 

Let us consider an autonomous Hamiltonian system with N degrees of freedom H(p, q) where p,q£ M N 
are the momenta and variables vectors. The Hamiltonian vector field may be written as 



x = jv x n = w( x ) , 



(i) 



where x 



G 



n2JV 



and 



J 



On —In 
Ijv On 



(2) 



is the standard symplectic matrix, being ljy the N x N identity matrix and Ojv the JVxJV matrix whose 
entries are all zero. 

Many chaos in dicators like the FLI |Froeschle et all Il997i|. th e MEGNO [Cincotta etaD . |2003| |. the 
SALI |Skokoj . l2001^ and more recently the GALI [Skokos et al.l . l2007l | require the time evolution of deviation 
vectors to be computed. These vectors, S(t), satisfy the variational equations given by 



S(t) = D x WS{t) = JVlH8{t) 



(3) 



where - DyW is the Jacobia n matrix of the vector field W and V^.^ is the Hessian matrix of H. 

In [Libert et all l2010f ]. the GSI has been introduced to simultaneously integrate both systems of 
equations ([!]) and Q by means of a symplectic integrator. This method assumes that % may be split into 
two separately integrable parts. For instance when 



the variational equations 

Sr 



can be written as 

'5 



-VKB 








(4) 



(5) 



while the associated Hamiltonian K, is expressed as 

/C( P , q, S p , <5 q ) = \slV 2 ^A8 p + i$Jv2 a 5«5 q = 
The above introduced notations will be use throughout the paper. 



^(p,<5 p ) + fi(q,5 q ) 



(6) 
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2.2. Integrator 

In the present wor k, two symplecti c integrators have been tested. Both are based on the Campbell-Baker- 



np_L 

Hausdorff (CBH) [Bourbaki . fl~972l ] formula which ensures that one can find a general explicit integrator 



with n steps of the form 

S n ( T ) = e c ^ L A e dirL B . . . e c n rL Ae d n rL B = ^Lg ^ ^ 

whose coefficients q and di have to be carefully chosen to get the required precision, i.e. reach the integrator 
order m. Integrating % at order m means that we exactly evaluate e TLg where 

g = A + B + 0{T m ). (8) 

In Laskar et all |200~D ] . four classes of symmetric symplectic integrators have been presented: 



SABA2 n (T) = g c i T ^A e diTL B e c n TLA e d n TL B e c n +irLA e d n TLB qC„tLa ^L\tL b ^c\tL a 
S ABA2 n +l (t) = e c i TL A e diTL B ^ _ e c n+ iTLA e d n +irL B e c n+ irLA _ gdiT-Lg ^cytLa 

SBAB2n(T") = e dirLB e c 2 TL A e d 2TL B e c n+1 TL Ae d n +iTL B e c n+1 TL A e d 2 rL B e c 2 TLA e diTL B 

SBAB2n+i('7") = e dirLB e C2rLA ...e dn+lTLs e Cn+2TLA e dn+1TLB ...e C2TljA e dlTLB . 



(9) 



If e := is small enough, it is shown that one can find specific coefficients such that 

G = A + B + 0(r 2n e + T 2 e 2 ). (10) 

Obviously, since the approximation error depends on the weight e of the perturbation, these classes of 
integrators are not well suited for Hamiltonian systems that are not perturbations of integrable ones. 

In general, for the Hamiltonian /C a ssociated to variatio nal equations, the ratio |i3|/|*4| is not necessarily 



small. Hence, the method presented in Laskar et al.1 . 120011 ] is not suitable. To tackle this p roblem, we hav e 



chosen to use another class of symmetric and explicit symplectic integrators presented in [Yoshida . 1990] 



The latter does not take into account the importance of the perturbation and can be defined starting from 
the basic bloc given by the second order Stdrmer- Verlet/Leap Frog scheme, T2 n d(T) = e^ TLA e TLB e^ TLjv , as 
follows: 

r 4t h(r) = T 2nd ( 2 _ 1 2 i/3 T ) T 2nd (- 2 ^ 2l/3 T ) T 2nd ( 2 _ 2l/3 r ) 1 1 1 » 



T 6th (r) = T 4th ( 2 - J 2 ^rJ T 4th ^--A-^rj T 4th (^iTr) ■ ^ 

In this case, the error depends only on the time step r, as in Eq. 0. While this method turns out to 
be very efficient for the variational Hamiltonian IC, it does not take advantage of the structure of % as a 
perturbation of an integrable system. 



3. MEGNO 
3.1. Definition 

According to (Cincotta et al. , 2003| ] , the Mean Exponential Growth factor of Nearby Orbits is defined as 



5(s 



-sds 



(13) 



where 5(s) denotes the Euclidean norm of S(s). A more useful and stable indicator is given by the mean 
MEGNO, namely the time-average: 



(14) 



W hile Y(t) may not converge nor admit a limit for t — > oo, it has been proven by Cincotta et all 
2003] that the asymptotic value of Y provides a good characterization of the regular or chaotic nature of 
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orbits. Basically, lim^oo Y(t) = 2 for quasi-periodic orbits on an irrational torus for a non-isochronous 
system and for orbits close to stable periodic ones. In the limit case where the orbit coincides with a stable 
periodic orbit, Y(t) asymptotically reaches zero. Considering irregular orbits, Y(t) increases linearly with 
time, being the slope half of the first Lyapunov exponent. 



3.2. Computation 



The computation of the MEGNO and its time-average requires both integrals ([13]) and (|Mj) to be solved 
accurately. Different methods are available. 

A straightforward approach is based on the introduction of two auxiliary functions v and w such that 



v (t) = tY(t) and w(t) = tY(t) , 
whose time evolution is directly given by the following differential equations: 



»(t) 



2 S -lt 



and w{t) = Y(t) 



v(t) 



(15) 



(16) 



's(ty ~ 5 2 t 

Obvi o usly, fl3l) and (fT6l) have to be comp uted with the same integrator (see e.g. [Gozdziewski et all . 



l200l| . |Valk et al. . l2009l |. Hinse et al. . 2010( |). In this case, the time step used to integrate Eq. ([3]) is fixed 



by the integration of Eq. (I16p . However, the use of auxiliary functions is less efficient within a symplectic 
integration scheme (in fact, (I16p are not generally Hamiltonian equations). Indeed, it is well-known that 
symplectic integrators show very good energy conservation properties, allowing us to consider larger step 
sizes and limit energy loss. In light of these considerations, other alternatives have been studied. 

In p articular, using the d efinition of the MEGNO for discrete time d ynamical systems (i.e. maps) 
given by [Cincotta et al. . 2003 ]. it has been proposed in Breiter et al. . 2005 ] to compute Y(t) and Y(t) by 
observing that a fixed step size integrator can be considered as equivalent to a discrete time map. Hence 

S(t + h) 



Y 



[Breiter et al., 200 



,](* + /0 



Y 



[Breiter et al., 200. 



6] (t + h) 



t + h 
tY 



Br. 



"^Rrcitcr ct al., 2005] (^) ""T" 2 III ■ 
tcr et al., 2005] (0 ~t~ ^ ^[Breiter ct 



S(t) 

iL, 2005] (t 



and 



h) 



t + h 



(17) 



h being the integration step size. Let us note that these formulas can be obtained by using a simple 
rectangular quadratur e method to solve b oth integrals (|13p and (|14p . On the other hand, a mixed scheme 
has been proposed in Gozdziewski 120031 ] . that relies on the computation of the MEGNO by using the 
so-called trapezoidal rule to compute the integral in the definition (|13p and then the discrete time 
approximation for the mean MEGNO. 

In the present work, we develop further this idea: we propose to use the trapezoidal rule to compute 
both integrals (|13p and (|14p . First, we rewrite (|13|) as 

Y(t) = 21og<5(i) - - / log 5{s)ds, (18) 



then using the trapezoidal rule we get 

Y(t + h) = 



t + h 



Y{t) + 2t ± h ln S±^ + 0{h3) 



t + h 



5(t) 



and 



Y(t + h) 



t + h 



t+h 



Y(s)ds 



t + h 



[tY(t) + 0.5h(Y(t) + Y(t + h))] + 0(h 3 



(19) 



(20) 



1 Given a real- valued function / and an interval [a, b], one has [Press et all |2007| 

f(x)dx = 0.5(6 - a)[f(a) + /(&)] + 0((b - aff") . 



the second derivative being estimated on the interval [a, b]. 
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Fig. 1. Arnold web. Two-dimensional phase plane (Ii,l2) represented using mean MEGNO values at t = 10 7 time units 
(the values greater than three have been fixed to three). A set of 600 x 600 uniformly distributed initial conditions has been 
integrated with (with time steps equal to 0.01). Other initial conditions are fixed to Is = 1, (pi = <f>2 = <j>3 = and 

v — 0.007. The dashed line represents the I\ = 2I2 resonance. The box encloses a part of it and is analyzed in more detail in 
Fig. El 



Let us observe that the above formulas (|19h and (I20p improve the aforementioned ones corresponding to 
lower order approximations of the integrals defining MEGNO and mean MEGNO. 

In the following, Eq. (I19p and (I20p will be used to compute the MEGNO whenever we use a symplectic 
scheme, whereas Eq. (I16p will be used with RK4. Let us note that in the rest of the paper we will be 
interested in asymptotic values of mean MEGNO i.e. mean MEGNO values at the end of the integration 
process. In order not to introduce any bias in the computation of the MEGNO and mean MEGNO, 
initial deviation vectors 6(0) will always be randomly chosen with uniform probability in the appropriate 
hyper sphere. 



4. Arnold Diffusion 
4.1. Model 



As shown in Libert et all [20101], the GSI proves to be more efficient than non-symplectic schemes to 
correctly identify the behavior of a given orbit especially on dynamics acting on long time scales. In this 
section, we will show that this is particularly relevant in case of slow chaotic diffusion. To that purpose, 
we decided to consider the following Hamiltonian system described in [Lega et"aL 20031 ]: 



^Arnold (h, h, h 



+ h + v 



cos((pi) + cos((j)2) + cos(4>3) + 4' 



(21) 
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0.29 0.295 0.3 0.305 0.31 0.315 0.32 0.325 0.33 



Fig. 2. Smaller portion of Arnold web. Two-dimensional phase plane (71, ^2) represented using mean MEGNO values at 
t = 10 time units (the values greater than three have been fixed to three). A set of 600 x 600 uniformly distributed initial 
conditions have been integrated with 24th (with time steps equal to 0.01). Other initial conditions are fixed to 73 = 1 and 
4>l — $2 = 4>Z — 0- The parallelogram represents the region in which 100 initial conditions have been considered to compare 
the GSI to a non symplectic integration scheme. 



where actions I±, I2, 13 £ K and angles , (j>2 , <p3 € T are canonically conjugate variables and v is assumed 
to be a small parameter. 

Given the structure of (|2Tj) . $1 = Ii, $2 = h and ^3 = 1. Hence, each straight line 

hh + fah + h = 0, {k h k 2 , fa) £ Z 3 \ {0} 

on the two-dimensional plane (Ji, I2) represents a resonance. 

As illustrated in Fig. [Tj most relevant resonances are clearly visible on the plane (Jj, J2), to form 
the so-called Arnold Web. Mean MEGNO values are shown for a grid of 600 x 600 equally spaced initial 
conditions in this plane. Other initial conditions are ^3 = 1 and 4>i = 4>2 = 4*3 = 0, and the parameter v has 
been fixed to 0.007, as in the rest of this work Thi s value needs to be small in order to avoid resonances 
overlap. Besides, as pointed out in [Lega et alJ . l200^ . the smaller the perturbation, the slower the diffusion. 



The GSI has been used with T^h integrator with a fixed time step r = 0.01 and an integration time of 
10 7 . The dashed line highlights the h = 2 Jo resonance. This web of resonances is of particular interest. 
Indeed, it has bee n proven (see [Arnold! . 19631 ]) that Arnold diffusion exists along resonances. Moreover, in 



Lega et all [2003], a numerical proof of diffusion along resonances for this model is given. 



The analysis presented in Section 14.21 will be performed in the region delimited by 0.29 < I\ < 0.33 
and 0.14 < I2 < 0.18, centered on the I\ = 2I2 resonance. This small region is enclosed in the box shown 
in Fig. [TJ An enlargement of this box is presented in Fig. [2j Again, a grid of 600 x 600 equally spaced initial 
conditions has been numerically integrated using the GSI with T^h integrator up to 10 7 time units. Other 
initial conditions and parameters are the same as the ones used to produce Fig. [H In the following analysis, 
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we will consider several orbits around the top hyperbolic border (in brown in Fig. [2]) of the resonance where 
diffusion is actually confined. 

4.2. Analysis 

In this section, we will compare the results on the correct determination of regular or chaotic orbits behavior 
obtained with the GSI and a non-symplectic scheme. Through an analysis of the maximum relative errors 
on the energy, percentage of correctly identified orbits and CPU time, we will show that our symplectic 
scheme outperforms the non-symplectic one in the determination of the behavior of orbits. 

As it is necessary to consider long integration times and a lot of different initial conditions, it has 
been decided to use and compare fourth-order integrators. In fact, this turns out to be a good compromise 
between reliability of the numerical results, as measured in term of relative energy loss, and number of eval- 
uations of the vector fields, translated easily into required CPU time. On the one hand, we considered the 
well-known non-symplectic fourth-order Runge Kutta (RK4) integrator. While it is simple to implement, 
it is also robust and efficient. On the other hand, both symplectic integrators presented in Section \2. 21 have 
been tested. As far as only the Hamiltonian system f|21 j) is concerned, SABA and SBAB integrator classes 
outperform Yoshida integrator. Indeed, given that v is small, this system can be seen as a perturbation of 
an integrable system. Hence, the error (jlpp is smaller than Yoshida's one ©. However, SABA and SBAB 
classes performances become rapidly poor for the Hamiltonian system K. Arnold associated to the variational 
equations. It is mainly due to the structure of /C Arnold- As time increases, the weight of the B part, say the 
perturbation, of the Hamiltonian /CAmoid may become larger than the A part, or vice versa without any 
possible control, see Fig. This implies that the error on the energy increases too. Fortunately, Yoshida 
integrators do not depend on this kind of consideration and are then more suited for this application. As 
already stated in Section 12.21 Yoshida integrators do not take advantage of the structure of "^Arnold but 
it is conversely useful for the integration of variational equations. Hence, for the purpose of the present 
study, it has been decided to compare RK4 integrator to T2 n d and T^h- 

The comparison has been performed considering 100 orbits whose initial conditions (/i(0), ^(O)) are 
uniformly distributed around the top hyperbolic border of the resonance (the parallelogram in Fig. [2]). 
First, for each one of these orbits, a reference value of the mean MEGNO has been computed with T^h, 
a time step equal to r = 0.01 and an integration time of 2.5 x 10 6 . Afterwards, for time steps between 
0.01 and 1, the same orbits have been numerically integrated with both methods (the GSI and the non- 
symplectic scheme) in order to obtain the corresponding MEGNO values. Eventually, a comparison with 
the reference MEGNO values has been done, enabling us to classify orbits as correctly identified or not. 
At the same time, CPU times and maximum relative errors in energy have been stored. In order to reduce 




Fig. 3. Relative weight of A and B parts (respectively A and B), for both Hamiltonian systems "^Arnold an d ^Arnold- One 
generic orbit and the associated variational equations have been integrated with with a time step r = 0.01. The horizontal 
line corresponds to v — 0.007. 



8 Ch. HUBAUX et al. 



numerical truncation errors, it has been decided to achieve this study using quadruple precision. 

A relevant indicator to test the goodness of the numerical integration, is to compare the maximum 
relative errors in energy 

AE/E= max \Ei(t) - Ei(0)\/\Ei(0)\ , (22) 

0<i<2.5xl0 6 

i=l,. ..,100 

Ei(t) being the energy, i.e. the value of the Hamilton function, at time t on the i— th orbit. 

A second important indicator, related to the speed of the integration algorithm, is the CPU time 

100 

Tcpv = ^2 T cpu > ( 23 ) 

i=l 

Tq PU being the CPU time needed to integrate the i— th orbit and the associated deviation vector on the 
defined time span. 

Both indicators are reported in Fig. |4] as functions of the time step for the different integrators. It 
appears that T^h shows always smaller energy loss than RK4 integrator. Moreover, as time step increases, 
this difference becomes larger too. Also note that the maximum error becomes larger with RK4 than with 
?2nd beyond t ~ 0.25. That means that, even if T2 n d is only a second order integrator, it is more reliable 
than RK4 when using big time steps. Another advantage of the GSI is the relatively low required CPU 
time (Tcpu) in comparison to RK4. This is particularly important as we are considering lots of different 
initial conditions and long integration times. Obviously, the lower-order T2 n d asks less CPU time than T^. 

Moreover, the GSI correctly identifies more orbits than RK4 as time steps increase (see Fig. [5]). Indeed, 
mean MEGNO values computed by means of RK4 are wrong for regular orbits when the time step is greater 
than 0.1. The percentage of well identified regular orbits even reaches zero while, at the same time, T^h 
is still beyond 50%. This difference is less discernible for chaotic orbits, since a small drift from the orbit 
and/or tangent direction does not lead to completely different behaviors. However, in the following, we will 
use the total percentage, p, that presents a summary of both results. Eventually, let us point out that the 
lower-order T2 n( j integrator manages to identify correctly approximately the same percentage of orbits as 




Fig. 4. Maximum relative errors in energy AE/E (left panel) and CPU time Tqp\j (right panel) as a function of the time 
step, in logarithmic scale. The integration time for this analysis has been set to 2.5 x 10 6 time units. The comparison involved 
100 orbits whose initial conditions (/^(O), /2(0)) have been taken around the top hyperbolic border of the Jj = 2I2 resonance 
(see Fig. [2}. Symbols are : (red) * RK4 integrator, (blue) O the 4th order Yoshida integrator and (black) □ the 2nd order 
Yoshida integrator. Straight lines denote best linear fits (left panel). 
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Fig. 5. Percentages of correctly identified orbits with respect to time step in logarithmic scale. The comparison involve fOO 
orbits whose initial conditions (ii (0), ^2(0)) have been taken around the top hyperbolic border of the I\ = 2/2 resonance (see 
Fig- EJ. Solid lines and dashed lines represent respectively the identification of regular orbits and chaotic orbits. 



10 
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10 ' 



10 
x 
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Fig. 6. Efficiency index <f> with respect to time step. The comparison involved fOO orbits whose initial conditions (7i(0), ^2(0)) 
have been taken around the top hyperbolic border of the I\ = 2/2 resonance. 



Eventual l y, all these observations can be recombined into a single efficiency index. As introduced in 
(Libert et alll20ld |. 



:p|log 10 (A£/£)||log 10 (T C Pu) 



(24) 



enables us to quantify the efficiency of each method. The larger </>, the better the method. On the one hand, 
the percentage of correctly identified orbits must be as big as possible. On the other hand, relative error in 
energy (g [0, 1]) and CPU time must remain low. The evolution of <j> with respect to time step is shown in 
Fig. [6l It turns out that this index presents larger values for computations realized with the GSI, coupled 
to T^h- It results obviously from previous considerations. For relatively small times steps (r < 0.02), the 
non-symplectic scheme (RK4) behaves similarly to the GSI. After that, this method is quickly penalized 
by its energy loss and lower percentage of well identified orbits. Eventually, efficiency index for the GSI 
used with T2 n( j presents similar behavior to the ones of but on a lower level. It comes directly from its 
larger maximum relative error in energy. However, at time step r ~ 0.2, the GSI with T2 n d becomes more 
efficient than the non-symplectic method, due to the joint effect of better energy conservation and number 
of correctly estimated orbits. 

Let us observe that the same analysis has been performed in double precision. Results concerning the 
percentage of correctly identified orbits are exactly the same. Obviously, CPU times are proportionally 
smaller than the present ones. The main difference occurs with the evolution of maximum relative energy 
errors with respect to time step which are not perfect straight lines anymore. In particular, numerical 
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truncation errors appear for too small relative errors in energy. This leads to an inaccurate efficiency index. 
5. Conclusion 

In this work, we have shown that the GSI is a powerful tool to characterize the regular or chaotic behavior 
of orbits in Hamiltonian systems. It proves to be especially efficient for dynamics acting on long time scales 
like the Arnold diffusion problem analyzed here, as energy conservation properties of symplectic integration 
schemes are fully pointed up. 

As it only asks for the knowledge of a single deviation vector, the MEGNO indicator is less time 
consuming than other chaos indicators and should be used jointly with the GSI. To numerically compute 
both MEGNO and mean MEGNO, we have introduced a quadrature method based on the trapezoidal 
rule. 

Furthermore, it has been shown that a particular attention has to be paid to the choice of the symplectic 
integrator. One has to bear in mind that both Hamiltonian systems associated to the equations of motion 
and the variational equations can present very different structures. For this reason, we argue for the use of 
Yoshida's symplectic integrator. 

Finally, our analysis shows that the GSI outperforms non-symplectic schemes. Indeed, large time steps 
are allowed, smaller computation times are needed and the energy loss is more limited. For all these 
considerations, we claim that the GSI method with the MEGNO indicator and Yoshida integrator turns 
out to be a reliable and time sparing method to correctly determine the behavior of orbits associated to 
Hamiltonian systems. 
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